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1 The Symbolic Computation of Operators 

After defining the twr expansions 

N 

exp(A) = ]T A '/ il 
:=0 
N 

loeJ' ■ 1) = ^(-l)*' +1 A7i, 

t=i 

a computer algebra vf+em such as MACSYMA or Maple will quickly com- 
pute 

xp(A) exp (B)) = A + B + 0 (N + 1). 
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Here 0 (N + 1) denotes terms containing a product of N -f 1 or more A’s 
and/or B' s. This computation depends crucially upon the fact that AB = 
BA; for objects for which this not true, certain correction terms enter. For 
example, if A and B are matrices, then in general AB ^ BA and 

log(exp( A) exp (J8)) = A + B-±AB + ±BA + . . . . 

The symbolic computation of operator expansions » uoh - a c these differ s 
in a number of ways from the symbolic computation of expressions in com- 
muting variables. The papers in this volume consider various aspects of such 
computations. In this introduction, we first discuss some of the capabilities 
that prove useful when performing computer algebra computations involving 
operators. These capabilities may be broadly divided into three areas: the 
algebraic manipulation of expressions from the algebra generated by opera- 
tors; the algebraic manipulation of the actions of the operators upon other 
mathematical objects; and the development of appropriate normal forms 
and simplification algorithms for operators and their actions. We then con- 
clude the introduction by giving a little background and a brief description 
of the problems considered by the authors in this collection. 

1.1 Algebra of Operators 

Let E u E 2 , . . . , Em denote operators and a a number. Then E\ + E<i, E 2 E 1 , 
and aE\ are all operators of the same type. That is, a finite set of operators 
naturally generate an algebra of operators under addition and composition. 
Let R{25i, . . . , Em} denote this algebra. This is just the free associative 
algebra over R generated by the symbols E \ 9 . . . , Em- The first capability 
of a computer algebra system for operators, then, is that it support the 
algebraic operations of addition and composition of operators. 

The first issue this raises is how to represent operators and operations 
on them in a context which has already usurped most of the notation for 
an algebra of expressions. Is it possible to use some multiplication operator 
(typically “*”) for operators, or should one use another notation? Maple[2] 
uses and Macsyma[4] uses for composition, but juxtaposition for 
application. Mathematical] takes no special notice of this, but allows the 
use of juxtaposition for multiplication. (It thereby maps f (a) into the same 
expression as a*f; if you intend to convey the notion /(a), you type f [a].) 
In fact, the notation and semantics of operators has at best been a patch 
added on to conventional general-purpose algebra systems; by contrast, the 
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most effective computer implementation of operator algebras has been quite 
application specific, as for example, in the case of manipulation of pertur- 
bation series. 

Yet the need for combining operator notation with the generality of 
widely available computer algebra systems dictates that we seriously ad- 
dress the representation and manipulation of operators. A good test for 
any system is to take some simple and familiar concepts and see if they can 
be naturally expressed and manipulated. For example, a operator that ex- 
presses “differentiation with respect to x” should be natural. Extensions to 
represent the 2nd or the nth derivative should be natural. Mixed derivatives 
(with respect to x and y) should be possible: a natural operation might be 
the interchange of orders of differentiation, and the combination of common 
variable “powers”. Evaluation of the derivative at a point, say the origin 
(the un-obvious /'( 0)) is a challenge [2]. Because algebra systems are gen- 
erally built to be extensible by procedure definition, data-type extension, 
transformation-rule (or pattern-match) extension, and aggregations of data, 
it is natural for one to hope that new operators can be defined, and their 
significant properties encoded, in a natural form. 

There is a tension between expressiveness and precision. In different con- 
texts, we have different expectations. Should the algebra system be expected 
to treat identical notations in different ways? Consider the notation ( L + 
a)(y). This might mean (say if L is a derivative operator) dy/dx + ay. Yet in 
other circumstances we might prefer that a constant a be an operator equiv- 
alent to a function which returns a, regardless of its argument. In that case, 
( L -f a)(y) = Ly + a. If an algebra system requires unambiguous notation, 
it may be unnatural in nearly all mathematical contexts. For the two vari- 
ations above, Maple would use <L®x+a*x|x>@y and <L®x+a|x>Qy, respec- 
tively while Macsyma would require a declaration declare(a,opscalar) 
then (L+a)y for the first and would use (L+constop(a))y for the latter. A 
proposal and implementation for Macsyma by T. H. Einwohner [3] would 
use the notation L+a*I I y to mean L(y) + a*y (Unfortunately, the use of 
“I” is at odds with the Maple convention). Einwohner’s design is somewhat 
more regular in reflecting the duality between expressions involving opera- 
tors (which can be quite arbitrary), and the results of applying the operator 
expression to arguments. The use of parameters is an important component 
of the design. For example, f (a) |y is an alternative notation for f (y,a). 
Finally, he does not further overload the use of as non-commutative 
multiplication by using it for operator composition. 

Another issue is how to best structure the data and the manipulative 
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algorithms for handling expressions from free associative algebras. These 
issues have been examined since the earliest days of computer algebra sys- 
tems, and are perhaps the aspects of the computer algebra of operators 
which can be attacked entirely by conventional means; this usually consti- 
tutes a mapping into algebraic tree representations, where most algorithms 
can be treated as tree-traversals (as described in any text on data struc- 
tures). On the other hand, truly efficient manipulation may dictate more 
compact representations, including so-called string-polynomial manipula- 
tions, matrix encodings, or other schemes that have been devised for hand 
manipulation. 

1.2 Actions of Operators 

The usefulness of operators to scientists and engineers derives not from their 
properties as abstract algebraic objects but rather from their interesting 
actions on other mathematical objects. For example, matrices act as linear 
transformations on vector spaces, vector fields act as derivations on spaces 
of functions, the physicist’s raising and lowering operators act on the space 
of polynomials, the algebra generated by shift operators acts on the space 
of sequences, and the algebra of maps acts on the space of domains. 

This brings us to the second major requirement on a computer algebra 
system for operators. They must support the action of operators on objects 
from the natural domains of definition of the operators (and presumably the 
natural domains of the algebra system). For example, having formed the 
matrix expression E = A + A 2 / 2! + . . . + A 4 / 4!, it is useful to be able to 
apply E to vectors. The merging of the two notations leads to a complexity 
and variety of actions that is probably the single most important barrier to 
satisfactory implementations of operators in computer algebra systems. The 
operator E above certainly looks like a polynomial in A; for some purposes 
it can best be treated as a polynomial; in other contexts as in section 1, it 
certainly is not equivalent. 

There is a more subtle, yet related issue. There is no natural mathemat- 
ical notation to specify the various interlocking operations possible on the 
mixed domain of operators and operands. For example, consider an operator 
F, and its inverse which we will conventionally write as F^ 1 (even though 
it may have rather little to do with 1/F) and an expression F F~ x x. It is 
natural for us to apply simplification and “cancel” the operators, yielding 
x . It is not so obvious to us, or to the computer, though, whether, in more 
complex situations to apply or simplify. Is FFx computed faster if F is 
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“squared” first? Or„is F(Fx), where the parentheses suggest some ordering, 
preferable? 

Sometimes partial simplification is useful. Consider a definition of a new 
operator F :=< 21"= 1 x \ x, n > where we use Maple’s notation to indicate 
that x and n are the formal names used for two actual “parameters” to 
the operator. If G is another operator (presumably operating on positive 
integers) then F(G,3) is G(l) + G{ 2) + <5(3). Consider < F(I,k ) | k >, 
where / is an identity operator. At what time (if ever) would a computer 
algebra system realize that this is equivalent to < k(k + l)/2 | k >? 

What kind of syntax do we supply for the user to define such a simplifi- 
cation? Furthermore, how do we deal with an action of a particular operator 
on a domain which is not supported? For example, should the syntax de- 
scribing the ’ : 'm of shift operators acting on sequences be the same as the 
syntax describing the action of vector fields as derivations on the space of 
functions? 

How can the implementation of operators as algebraic objects be best 
merged with the implementation of operators as objects acting on other 
objects in the computer algebra system? 

It appears that linguistically, two approaches are possible, and these 
are not addressed in Maple or Macsyma. One is to require unambiguous 
specification of operations such as operator simplification (so they would 
occur on command), a distinct notation for application of operators, and an 
explicit translation (back and forth) from operator to functional notation. 

Another approach is to use the domains of objects (themselves possi- 
bly other operators) to disambiguate the meanings of operators, at least 
whenever possible. This may require fairly sophisticated pattern-matching 
support which checks the types of arguments. An example used by Gonnet 
[2] illustrates this. Consider the expression ax D x D x y. If we assume D is 
an operator, each of the three “multiplications” may be a different concept. 
The first is symbolic multiplication by a constant. The second is composi- 
tion, and the third is application. Yet we were able to disambiguate this by 
looking at the types of the “arguments” to x. A constant on the left, a, 
means traditional multiplication; a non-operator on the right, y means that 
the operator to the left is being applied to it. A multiplication between two 
operators means composition. Note that it would be an error to simply work 
from the right to the left, applying as we go, although for this expression it 
might work. Consider a non-integrable form z, and the integration operator 
D~ l . Then D D~ l z could not be simplified, because the application of the 
integration operator would not “work” (unless D were specially coded to 
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“unwrap” an integral). 

There are many open questions, and it appears to us that the best pro- 
cess for making headway in the introduction of operator calculi in algebraic 
manipulation is to explore serious applications, ones which are challenges for 
humans and computers. Without further experience such as we see in this 
volume of papers, it is too easy to make false starts and empty promises. 
We do not know, for example, if it make sense for the language describing 
the action of matrices on vectors to be the same as the language describing 
the action of differential operators on the space of functions. 

While Maple proposes a methodology primarily based on an extension 
of procedures, with some original notation for operators, Macsyma uses (for 
the most part) already existing notation for non- commutative multiplica- 
tion; we expect that users of the Mathematica system will primarily use 
pattern matching and notations which look fundamentally like functional 
application. Each approach has its advantages and advocates. 

1.3 Normal Forms and Simplification 

Data structures and algorithms to manipulate operators depend crucially 
on the action of the operators upon objects from the natural domains of 
definition of the operators. Normal forms for expressions built from matrices 
are probably not the best normal forms for expressions built from maps. 
Questions about normal forms, simplification, and efficient evaluation of 
operator expressions are by and large open. 

It appears that another significant area of application of a computer 
algebra system is the manipulation of operator expressions to produce sim- 
plified or if possible normal forms for the various types and combinations 
of supported operators by the computer algebra system. This may involve 
transformation of operator expressions to a well-supported domain (typically 
polynomials), or collections of rules which drive expressions algorithmically 
or heuristically, toward equivalent but simpler forms. 

Given the state of the art, it seems inevitable that a successful general- 
purpose system will have to provide some facility for users to implement 
their own normal forms and simplification algorithms for more specialized 
types of operators and actions. 
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Examples of Operators and their Domains 

In this section we give brief descriptions of the computer algebra computa- 
tions that arise when working with various operators and their actions. The 
papers in this collection give careful and complete descriptions of how the 
authors dealt with some of the issues mentioned above. 


2.1 Perturbation Theory 


A classical problem in perturbation theory is to compute approximate so- 
lutions to differential equations containing small parameters. Consider van 
der Pol’s equation 


d 2 x 

dF 


z o, dz 

+I _ f( i _ I )_ 


= 0 , 


where e ^ 0 is a small parameter. The starting point is to expand a putative 
solution t —* x(t) in a power series in e 


x(t) = x 0 (t) + exi(t) + e 2 x 2 (0 + • • • , 


and then substitute this series into the original differential equation to ob- 
tain a sequence of differential equations (one for each power of e) for the 
unknown coefficient functions x,(f). These auxiliary differential equations 
have the property that the ith equation involves only the coefficient func- 
tions xo(t), . . ., x,(t), so that the sequence of differential equations can be 
solved recursively. In the example above, the equations are 

d 2 

x 0 (t) + x 0 (t) = 0 

■p*i(i) + a:i(f) = (^*o(0) (! “ *o 2 (0) 


There are several ways of approaching these types of computations. Let- 
ting x(t) = y(t), Van der Pol’s equation can be written as the first order 
system 


*(0 = lK0 

y(t) = -x(t) + e(l - x 2 (t))y(t). 
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Let E\ denote the vector field yd/dx — xd/dy and let E 2 denote the vector 
field (1 — x 2 )yd/dy. Then Van der Pol’s equation becomes 

z(t) = (E \+eE 2 )(z(t)), 

where z(t) = (x(t), y(t)). Notice that E 2 E\ ^ E\E 2 . 

Let T„ denote the operator which acts upon functions t — ► z(t ) and 
returns the first nth terms in a power series expansion in e of the function. 
Then the ith auxiliary equation is equivalent to the system of equations 
which is the coefficient of e‘ in the expansion of 

(d/df - (Ei + eE 2 ))T n (z(t)) = 0, 

for n sufficiently large. From this point of view, perturbation theory is 
concerned with the symbolic algebra of noncommutative operators such as 
Ei and E 2 acting on the domain of power series expansions of the form 

f{x,y) = fo{x,y) + €/i(*,2/) + t 2 f2{x,y) + ■ • • • 

A different, but related point of view is used by R. Rand in his contribu- 
tion “Perturbation Methods and Computer Algebra”. In this paper Rand 
describes a computer algebra system built using MACSYMA, which in a sys- 
tematic fashion can perform serious perturbation computations, including 
the computation of normal forms and center manifold reductions. 

2.2 Finite Difference Operators and Domains 

Consider the heat equation in a bounded region Cl of the Euclidean plane 
du d 2 u d 2 u 

1 e - 

To compute the temperature u(x,t ) numerically using finite differences, we 
need to discretize the domain Cl, the function u(x,t), and the differential 
operator 

d d 2 d 2 

dt dx 2 dy 2 

This can be done in many ways. Let D denote a finite mesh of points 

xi = xq + iAx for i = 0, . . . , / 

Vj = yo+j&y for j = 0,...,J 
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covering the region, ft, and let u 7 }- denote the temperature at time nAt at 
the mesh point (zo + iAx,yo + jAy ). 

With the shift operators 

E *( k K,j = u i-k,j 

E y( k ) u i,j = 

we can define the difference operators 

= (^x(-l) - 2£*(0) + £?*(1)) 

= ^,.,--2 vtj+vf-u 
- -2E y (0) + E y (l))ul J 


and compute the temperature u"* 1 given the temperature ujj implicitly 
using the scheme 


= <j 


+ 


5 (&)(«< 


+ 1 

'* J 


+ 


+ 


s l<? 


+ *}<i) ■ 


Notice that the basic objects are: domains, such as D; functions defined 
on domains, sik aerators defined on the functions, such as 

id ft 2 tie latter two objects can both be thought 

> il Vj it u. domain can be thought of as a map from the do- 
main to the range of the function; an operator on functions can be thought 
of as a map from the space of functions on domains to itself. In the pa- 
per U FIDIL: A Language for Scientific Programming” by P. Hilfinger and 
P. Colella, the language FIDIL (Finite Difference Language) is introduced 
which makes domains and maps basic data types and provides for the effi- 
cient computation of objects built from them. This makes the translation 
of standard numerical algorithms into programming statments very simple. 

Related work is contained in [lb 


2.3 Automatic Der : of Dynamical Equations 


The time evolution 
bodies joined tor 
complicated. T 
can be difficm 


anical system consisting of a system of rigid 
massless ball and socket joints can be quite 
to write down the correct equations of motion 
. ^uld be very useful if a program could automatically 
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derive the equations of motion of this type of mechanical system. Similarly, 
an electrical circuit consisting of resistors, capacitors, and voltage sources 
can exhibit interesting dynamical behavior. It is an interesting problem 
to write a program whose input consists of a description of a mechanical 
or electrial system and whose output consists of the differential equations 
governing the time evolution of the state variables of the system. 

A description of the system would include the following: 

System parameters . System parameters must be defined and specified. For 
example, the moment of inertia of a body is defined to be f body Q • Qdm(Q), 
where m(Q) is the mass distribution of the body. To describe the system, 
the mass distributions and topology of the connections of the rigid bodies 
must be given, and the moments of inertia must be computed. 

Dynamical variables . The dynamical variables must be defined. For exam- 
ple, the rotation matrix B(t ) of a rigid body, which specifies how the body 
is turning in space, must be defined. Using this, the angular velocity £l(t) 
of the body can be defined via il(t) = B(t)B~ l (t). As a another example, 
Newton’s Law F(t) = dp(t)/dt defines the force F(t ) acting on a body in 
terms of the momenta p(t) of the body. Both these example are typical, in 
the sense that dynamical variables are often defined by differentiating other 
dynamical variables. Notice that this gives rise to differential identities sat- 
isfied by the dynamical variables. 

Algebraic relations. The dynamical variables not only satisfy differential 
relations, but typically algebraic relations determined by the geometry of 
the particular system. For example, if ri(t) and r 2 (t) denote the positions 
of the center of mass of bodies 1 and 2 respectively, r(t) denotes the position 
of the center of mass of the ensemble of the two bodies, di(£) and d 2 (t) their 
initial displacements, and B\(t) and B 2 {t) the rotation matrices describing 
the rotation of the bodies, then 

r(t) = Bi(t)di(t) + B 2 (t)d 2 (t). 

State variables. Because of the algebraic and differential relations satisfied 
by the dynamical variables, it is possible to select some of the dynamical 
variables, called state variables, and to construct all of the other dynamical 
variables from these. For example, for two bodies, only the position and 
velocity of the two bodies and their center of mass need be specified. Alter- 
nately, the position and momenta of the two bodies, and of their center of 
mass, can be specified. 
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Evolution of state variables. Finally, differential equations giving the time 
evolution of the state variables need to be given, together with the corre- 
sponding initial conditions. For example, the initial positions and velocities 
of each of the bodies and of each of the joints must be specified in order for 
the dynamical equations giving the time evolution of the state variables to 
be integrated. 

The contribution “Multibody Simulation in an Object Oriented Pro- 
gramming Environment” by N. Sreenath and P.S. Krishnaprasad undertakes 
the automatic derivation of the equations of motion for systems of coupled 
rigid bodies ir the plane, while the contribution “The Dynamicist’s Work- 
bench- T ‘ Reparation of Numerical Experiments” by H. Abelson 

an< electrical networks consisting of resistors, capacitors, 

clli* .) ouu*tCS. Both papers start with a description of the system, 

and oy using a variety of symbolic and symbolic/numeric techniques even- 
tually compute numerical simulations of the differential equations which the 
programs automatically compute. 

The papers treat the problem of finding state variables differently. The 
paper by Abelson and Sussman finds state variables using a symbolic Newton 
iteration to eliminate relations among dynamical variables. The paper by 
Sreenath and Krishnapr: Mses symmetries of the mechanical system and 

a technique from geomeu chanics called reduction to derive equations 
of motion on a smaller dimensional phase space consisting of state variables. 

Studies such as these make heavy use of symbolic computation of oper- 
ators. For example, the dynamical equations of a system of coupled rigid 
bodies can be written as F(t ) = {F, H} y where H is the Hamiltonian of the 
system, and both F and H evolve on the appropriate phase space P. Here 
{•,•} is a noncommutative operator defined on 

C°°(P) x C°°(P) 

2.4 L ; ° ~ s and Vector Field Systems 

Recaj] , a smooth vector field F(x ) defined in a neighborhood of the 
origin in is of the form 

F(X) = £>(X)A 

i=l ax * 

where aj(x), . . . ,a^(x) are smooth functions defined in a neighborhood of 
the origin. The Lie brackets of two vector fields F(x) and G(x ) is defined 


C°°(P). 
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to be 


[F,G] = F(x)G(x)-G(x)F(z). 

It is a fundamental property of Lie brackets that although they appear to 
be second order differential operators they are in fact first order differential 
operators because of the cancellation of the second order partial derivatives. 
In other words, the Lie bracket of two vector fields is a vector field. It 
therefore makes sense to form iterated Lie brackets such as [[F,G],G]. 

Observe that the Lie bracket of two vector fields is a measure of how far 
from commuting the two vector fields are. For example, the Lie brackets 
of the coodinate vector fields F{ = d/dx{ are all zero, indicating that these 
vector fields all commute with each other. In other words, F{Fj = FjFi. 
On the other hand, the Lie bracket of the vector fields F\ = x 2 d/dx\ and 
F 2 = d/dx 2 is equal to 

[F 2 ,F 1 ]^d/dx l . 

Let F(x) and Gi(x), . . . , G m (x) denote smooth vector fields defined in 
a neighborhood of the origin in R^. It turns out that the local behavior of 
the nonlinear control system 

m 

x(t) = F(x(t )) + ^2 Ui(t)Gi(x(t)), x(0) = x 0 6 R n , 

i=i 

is determined by the algebraic properties of iterated Lie brackets of the F 
and G’s. This is analogous to the fact that the local behavior of a smooth 
function is determined by its Taylor coefficients. Because of this fact, ques- 
tions about the dynamical behavior of control trajectories can be reduced 
to symbolic questions about the algebraic properties of the noncommutative 
operators F, G\ , . . . , G m acting on the domain C°°(R^) of smooth functions 
on R*. 

The contribution “Symbolic Computations in Differential Geometry Ap- 
plied to Nonlinar Control Systems” by 0. Akhrif and G. L. Blankenship 
describes a software package written in the computer algebra system MAC- 
SYMA, which through the symbolic computation of the proper Lie bracket 
expressions, can analyse system theoretic properties of a nonlinear control 
system, such as feedback equivalence, left-invertibility, or the design of con- 
trol laws. 

The contribution “Vector Fields and Nilpotent Lie Algebras” by M. 
Grayson and R. Grossman considers conditions on the Lie brackets which 
guarantee that the trajectories of a vector field system can be explicitly in- 
tegrated in closed form. In general this cannot be expected and the game 
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is to find as large a_class of such systems as possible. The paper gives some 
simple examples of such systems which were computed used the computer 
algebra system Maple. 
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